library(ggplot2)
library(foreign)
setwd("~/Replication Code/")

##---Loading in data ----
anes9297 <- read.dta("./data/anes9297_subset.dta")

ccap08 <- read.dta("./data/ccap08_subset.dta")

ccap12 <- read.dta("./data/ccap12_subset.dta")

vsg.dat <- read.dta("./data/vsg_subset.dta")

ccap16 <- read.dta("./data/ccap16_subset.dta")

##--Models----
#---ANES MODELS 92-94 -----
dat1 <- na.omit(anes9297[c("white_int92", "white_int94",
                            "pid92_7scaled", "pid94_7scaled",
                            "rr92_scaled", "rr94_scaled", "V940005")])

### PID on RR
m1 <- lm(rr94_scaled ~  pid92_7scaled  + rr92_scaled, 
         data = dat1, weights = V940005, 
         subset = white_int92 == white_int94)
summary(m1)

### RR on PID
m10 <- lm(pid94_7scaled ~  pid92_7scaled  + rr92_scaled, 
          data = dat1, weights = V940005, 
          subset = white_int92 == white_int94)
summary(m10)

#---CCAP 2008 -----
dat2 <- na.omit(ccap08[c("pid7r_m_sc", "rr_m_scaled", "rr_o_scaled", "pid7r_o_sc", "weight")]) 

### PID on RR
m2 <- lm(rr_o_scaled ~ pid7r_m_sc  + rr_m_scaled, data = dat2, weights = weight)
summary(m2)

### RR on PID
m20 <- lm(pid7r_o_sc ~ pid7r_m_sc  + rr_m_scaled, data = dat2, weights = weight)
summary(m20)

#---CCAP 2012: P2----
ccap.p2 <- subset(ccap12, p2 == 1, 
                  select=c("rr_p1_scaled", "rr_p2_scaled", "pp_pid7r_sc", "pid7_sc", "weight"))
dat3 <- na.omit(ccap.p2)
# PID on RR
m3 <- lm(rr_p2_scaled ~  pp_pid7r_sc  + rr_p1_scaled, 
         data = dat3, weights = weight)
summary(m3)

# RR on PID
m30 <- lm(pid7_sc ~  pp_pid7r_sc  + rr_p1_scaled, 
          data = dat3, weights = weight)
summary(m30)

#---CCAP 2012: P3-----
ccap.p3 <- subset(ccap12, p3 == 1, 
                  select=c("rr_p1_scaled",  "rr_p3_scaled",  "pp_pid7r_sc", "pid7_sc", "weight"))
dat4 <- na.omit(ccap.p3)
# PID on RR
m4 <- lm(rr_p3_scaled ~  pp_pid7r_sc  + rr_p1_scaled, data = dat4, weights = weight)
summary(m4)

# RR on PID
m40 <- lm(pid7_sc ~  pp_pid7r_sc  + rr_p1_scaled, 
          data = dat4, weights = weight)
summary(m40)

#---VSG 2012-2016 ----
dat5 <- na.omit(vsg.dat[c("rr_sc_12", "rr_sc_16", "pid7_sc_12", "pid7_sc_16", "weight")])

# PID on RR
m5 <- lm(rr_sc_16 ~ pid7_sc_12 + rr_sc_12, data = dat5, weights = weight)
summary(m5)

# RR on PID
m50 <- lm(pid7_sc_16 ~ pid7_sc_12 + rr_sc_12, data = dat5, weights = weight)
summary(m50)

#---CCAP 2016 -----
dat6 <- na.omit(ccap16[c("b_rr_sc", "p_rr_sc", "b_pid7_sc", "p_pid7_sc", "weight_post")])


# PID on RR
m6 <- lm(p_rr_sc ~ b_pid7_sc + b_rr_sc , 
         data = dat6, weights = weight_post)
summary(m6)

# RR on PID
m60 <- lm(p_pid7_sc ~ b_pid7_sc + b_rr_sc, 
          data = dat6, weights = weight_post)
summary(m60)

##---Helper function-----
to_plot <- function(model = NULL, vars = NULL, intercept = TRUE, robust = FALSE, ci = 95, shift = "range", num_sd = 1){
  ## Collect Variables to Plot
  if(length(vars) == 0){
    vars <- names(coef(model)) 
    vars <- grep("(Intercept)", vars, value = T, invert = T)
  }
  ## Get Coefficients
  if(intercept == TRUE){
    vars <- c("(Intercept)", vars)
    m_coefs <- coef(model)[vars] 
  }
  if(intercept == FALSE){
    m_coefs <- coef(model)[vars] 
  }
  ## Get Standard error
  if(robust == FALSE){
    m_ses <- sqrt(diag(vcov(model)))[vars] 
  }
  if(robust == TRUE){
    m_ses <- sqrt(diag(sandwich::vcovHC(model, type = "HC1")))[vars] 
  }
  
  ## Get CI
  # check ci format
  if(ci > 1){
    ci <- ci/100
  }
  alpha <- 1 - ci
  # store critical value
  ci_z <- abs(qnorm(alpha/2))
  ## Get range for variables in data
  tmp_dat <- model.matrix(model)
  
  shift <- tolower(shift)
  if(shift == "range"){
    if(length(vars) == 1){
      dat_shift <- diff(range(tmp_dat[,vars])) 
    }
    if(length(vars) > 1){
      dat_shift <- diff(apply(tmp_dat[,vars], 2, range)) 
    }
    if(intercept == TRUE){
      dat_shift[, "(Intercept)"] <- 1
    } 
  }
  ## Get SD for variables if this is desired
  if(shift == "sd"){
    if(length(vars) == 1){
      dat_shift <- sd(tmp_dat[,vars])*num_sd
    }
    if(length(vars) > 1){
      dat_shift <- apply(tmp_dat[,vars], 2, sd)*num_sd
    }
    if(intercept == TRUE){
      dat_shift["(Intercept)"] <- 1
    }
  }
  
  # Combine 
  m_coef_shift <- m_coefs*dat_shift
  m_ses_shift <- m_ses*dat_shift
  m_lci <- m_coef_shift - ci_z*m_ses_shift
  m_uci <- m_coef_shift + ci_z*m_ses_shift
  
  tmp_df <- data.frame(var = vars, coef = as.numeric(m_coef_shift), lci = as.numeric(m_lci), uci = as.numeric(m_uci))
  return(tmp_df)
}

##---Figure 1-------
CI <- 84
anes_rr <- to_plot(m1, vars = "pid92_7scaled", intercept = F, robust = T, shift = "sd", ci = CI)
anes_rr$df <- "ANES"
anes_rr$dv <- "rr"
ccap08_rr <- to_plot(m2, vars = "pid7r_m_sc", intercept = F, robust = T, shift = "sd", ci = CI)
ccap08_rr$df <- "ccap08"
ccap08_rr$dv <- "rr"
ccap12_p2_rr <- to_plot(m3, vars = "pp_pid7r_sc", intercept = F, robust = T, shift = "sd", ci = CI)
ccap12_p2_rr$df <- "ccap12_p2"
ccap12_p2_rr$dv <- "rr"
ccap12_p3_rr <- to_plot(m4, vars = "pp_pid7r_sc", intercept = F, robust = T, shift = "sd", ci = CI)
ccap12_p3_rr$df <- "ccap12_p3"
ccap12_p3_rr$dv <- "rr"
vsg_rr <- to_plot(m5, vars = "pid7_sc_12", intercept = F, robust = T, shift = "sd", ci = CI)
vsg_rr$df <- "vsg"
vsg_rr$dv <- "rr"
ccap16_rr <- to_plot(m6, vars = "b_pid7_sc", intercept = F, robust = T, shift = "sd", ci = CI)
ccap16_rr$df <- "ccap16"
ccap16_rr$dv <- "rr"

anes_pid <- to_plot(m10, vars = "rr92_scaled", intercept = F, robust = T, shift = "sd", ci = CI)
anes_pid$df <- "ANES"
anes_pid$dv <- "pid"
ccap08_pid <- to_plot(m20, vars = "rr_m_scaled", intercept = F, robust = T, shift = "sd", ci = CI)
ccap08_pid$df <- "ccap08"
ccap08_pid$dv <- "pid"
ccap12_p2_pid <- to_plot(m30, vars = "rr_p1_scaled", intercept = F, robust = T, shift = "sd", ci = CI)
ccap12_p2_pid$df <- "ccap12_p2"
ccap12_p2_pid$dv <- "pid"
ccap12_p3_pid <- to_plot(m40, vars = "rr_p1_scaled", intercept = F, robust = T, shift = "sd", ci = CI)
ccap12_p3_pid$df <- "ccap12_p3"
ccap12_p3_pid$dv <- "pid"
vsg_pid <- to_plot(m50, vars = "rr_sc_12", intercept = F, robust = T, shift = "sd", ci = CI)
vsg_pid$df <- "vsg"
vsg_pid$dv <- "pid"
ccap16_pid <- to_plot(m60, vars = "b_rr_sc", intercept = F, robust = T, shift = "sd", ci = CI)
ccap16_pid$df <- "ccap16"
ccap16_pid$dv <- "pid"

preds <- rbind(anes_rr,ccap08_rr, ccap12_p2_rr, ccap12_p3_rr, vsg_rr, ccap16_rr,
               anes_pid,ccap08_pid, ccap12_p2_pid, ccap12_p3_pid, vsg_pid, ccap16_pid)
preds$ci <- 84

CI <- 95
anes_rr <- to_plot(m1, vars = "pid92_7scaled", intercept = F, robust = T, shift = "sd", ci = CI)
anes_rr$df <- "ANES"
anes_rr$dv <- "rr"
ccap08_rr <- to_plot(m2, vars = "pid7r_m_sc", intercept = F, robust = T, shift = "sd", ci = CI)
ccap08_rr$df <- "ccap08"
ccap08_rr$dv <- "rr"
ccap12_p2_rr <- to_plot(m3, vars = "pp_pid7r_sc", intercept = F, robust = T, shift = "sd", ci = CI)
ccap12_p2_rr$df <- "ccap12_p2"
ccap12_p2_rr$dv <- "rr"
ccap12_p3_rr <- to_plot(m4, vars = "pp_pid7r_sc", intercept = F, robust = T, shift = "sd", ci = CI)
ccap12_p3_rr$df <- "ccap12_p3"
ccap12_p3_rr$dv <- "rr"
vsg_rr <- to_plot(m5, vars = "pid7_sc.12", intercept = F, robust = T, shift = "sd", ci = CI)
vsg_rr$df <- "vsg"
vsg_rr$dv <- "rr"
ccap16_rr <- to_plot(m6, vars = "b_pid7_sc", intercept = F, robust = T, shift = "sd", ci = CI)
ccap16_rr$df <- "ccap16"
ccap16_rr$dv <- "rr"

anes_pid <- to_plot(m10, vars = "rr92_scaled", intercept = F, robust = T, shift = "sd", ci = CI)
anes_pid$df <- "ANES"
anes_pid$dv <- "pid"
ccap08_pid <- to_plot(m20, vars = "rr_m_scaled", intercept = F, robust = T, shift = "sd", ci = CI)
ccap08_pid$df <- "ccap08"
ccap08_pid$dv <- "pid"
ccap12_p2_pid <- to_plot(m30, vars = "rr_p1_scaled", intercept = F, robust = T, shift = "sd", ci = CI)
ccap12_p2_pid$df <- "ccap12_p2"
ccap12_p2_pid$dv <- "pid"
ccap12_p3_pid <- to_plot(m40, vars = "rr_p1_scaled", intercept = F, robust = T, shift = "sd", ci = CI)
ccap12_p3_pid$df <- "ccap12_p3"
ccap12_p3_pid$dv <- "pid"
vsg_pid <- to_plot(m50, vars = "rr_sc_12", intercept = F, robust = T, shift = "sd", ci = CI)
vsg_pid$df <- "vsg"
vsg_pid$dv <- "pid"
ccap16_pid <- to_plot(m60, vars = "b_rr_sc", intercept = F, robust = T, shift = "sd", ci = CI)
ccap16_pid$df <- "ccap16"
ccap16_pid$dv <- "pid"

preds2 <- rbind(anes_rr,ccap08_rr, ccap12_p2_rr, ccap12_p3_rr, vsg_rr, ccap16_rr,
                anes_pid,ccap08_pid, ccap12_p2_pid, ccap12_p3_pid, vsg_pid, ccap16_pid)
preds2$ci <- 95

preds <- rbind(preds, preds2)

preds$df <- factor(preds$df,
                   levels = c("ANES", "ccap08", "ccap12_p2", "ccap12_p3", "vsg", "ccap16"),
                   labels = c("1992-1994\nANES", "2008\nCCAP", "2012\nCCAP\nMarch", 
                              "2012\nCCAP\nAugust", "2012-2016\nVOTER\nSurvey", "2016\nCCAP"))
preds$cl <- factor(preds$dv,
                   levels = c("rr", "pid"),
                   labels = c("Lagged Partisanship", "Lagged Racial Resentment"))
preds$dv <- factor(preds$dv,
                   levels = c("rr", "pid"),
                   labels = c("Predicted\nRacial Resentment", "Predicted\nPartisanship"))


ggplot(preds, 
       aes(x = df, y = coef, group = cl)) +
  geom_linerange(data = subset(preds, ci == 95),
                 aes(ymin = lci, ymax = uci), color = "grey", 
                 position = position_dodge(width = .5), lwd = .75) +
  geom_linerange(data = subset(preds, ci == 84),
                 aes(ymin = lci, ymax = uci), 
                 position = position_dodge(width = .5), lwd = 1.25) +
  geom_point(size = 3, aes(shape = cl), 
             position = position_dodge(width = .5)) +
  ggthemes::theme_hc() +
  scale_shape_manual(values = c(15, 16)) +
  scale_y_continuous(labels = numform::ff_num(zero = 0, digits = 2),
                     limits = c(-.01, .07)) +
  theme(legend.title = element_blank(),
        legend.text = element_text(size = 16),
        axis.text.x = element_text(size = 16),
        axis.text.y = element_text(size = 16),
        axis.title = element_text(size = 18, face = "bold"),
        strip.background = element_blank(),
        strip.text.x = element_text(size = 16),
        plot.title = element_text(size = 18, hjust = 0.5),
        plot.margin = unit(c(0, 0, 0, 0), "cm")) +
  labs(title = "",
       y = "Estimated Effect",
       x = "")
# ggsave("./figures/comb_sdshift.pdf", width = 8, height = 6)